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ABSTRACT 

Full NVT quantum statistics of the ion is simulated at low densities using the path integral 
Monte Carlo approach. For the first time, the molecular total energy, partition function, free energy, 
entropy and heat capacity are evaluated in temperatures relevant for planetary atmospheric physics. 
Temperature and density dependent dissociation recombination reaction balance of the molecule and 
its fragments above 4000 K is described, and also, the density dependence of thermal ionization above 
10000 K is demonstrated. We introduce a new well-behaving analytical model for the molecular parti- 
tion function of the H^^ ion for the temperature range below dissociation and fit the parameters to the 
energetics from our simulations. The approach presented here can be regarded as an extension of the 
traditional ab initio quantum chemistry beyond the Born-Oppenheimer approximation to description 
of nonadiabatic phenomena, and even further, account of nuclear quantum dynamics. 
Subject headings: molecular data - stars: atmospheres 



1. INTRODUCTION 



The molecular ion has been a subject of a num- 
ber of theoretical and expe rimental studie s since its 
first experimental detection (|Thomsonlll911| ). Because 
of its rapid formation through the exothermic reaction 
{AE w -1.7 eV) 



H2 + H+ ^ H+ + H 



(1) 



the H3 ion is expected in a ny active environment con - 
taining molecular hydrogen (jNeale and Tennvsonlll995l ). 
and thus , it is encou ntered, e.g. in the studies of the giant 
planets (jLvstrup et al. 2008; Koskinen et al. 2009). 

In planetary atmospheric physics, importance of the 
ion lies in its c apability to act as a cooling agent via 
infrared radiation (iNeale et al\\\M^ iHarris et a^] 120041 : 
iKoskinen et al\ 120071 ). The atmospheric models taking 
into account this cooling are commonly based on the high 
temperature molecular p artition function of the H;^ ion 
(jNeale and Tennvson|[l995i ). Evaluation of the partition 
function faces, however, a few challenges of which the 
first one is finding a good approximation to the infinite 
summation over all rovi brational quantum states w ith ac- 
curate enough energies (jNeale and Te nnvson 19951 This 
has usually been worked out with calculations of a finite 
number of stat es from, e.g. a sem i-empirical potential 
energy surface pinelh et ai.lll995[ ). 

The next challenge comes with the changing geometry 
of the ion at finite temperature. The rovibrational 
model needs to be extended for calculations of correct 
energetics for the emerging linear g eometry of the weakly 
bound molecule (jNeale aZ.lll996( ). 

Finally, at about the same temperature, where linear 
geometries start contributing, the molecule may also dis- 
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sociate to its fragments, and in fact, the equilibrium 



H+ 0H2 

o H+ 



P 
H 



o 2H + p ' 
o H + 2p+ + e" 
o 3p+ + 2e" 



(2) 



needs to be considered above about 4000 K, the balance 
depending strongly on both the temperature and the den- 
sity. ^ ^ 

This brings forth two questions, at the least. First, 
how relevant it is to consider the molecular energetics 
and related partition function at temperatures where the 
molecule has dissociated and appears in form of frag- 
ments of the equilibrium reaction, Eq. only. Sec- 
ondly, the balance of the equilibrium reaction may be 
strongly affected, not only by the density, but also by the 
environment including the neutralizing negative counter- 
parts of the positive Hj^. Thus, the thermal dissociation- 
recombination balance above about 4000 K gives rise to 
problems, which have not been taken into account in this 
context, yet. 

In this study, using the path integral quantum Monte 
Carlo (PIMC) method we have carried out the first sim- 
ulations of the full quantum statistics of the ion, 
described by Eq. ([2|), at low densities and high tem- 
peratures ranging from 160 K up to about 15000 K. 
PIMC is the method to meet the above challenges: we 
need not make any approximations or restrictions in the 
summing over states, geometries or quantum description 
of dynamics. The finite temperature is inherent in the 
PIMC approach and the Coulomb many-body treatment 
of the particle interactions is exact. The PIMC method 
is computational ly expensive, but feasible for smal l 
enou g h syst ems. ('Ceperlevlll995l: |Pierce and Manousaki; 
1999r iKwon and Whalcv" 1999: iKnoh and Marx 200" 
Cuervo and RoMi2006t iKvlanpaa and Rantalall2009i V 

The conventional quantum chemical ab initio descrip 
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tion of the ion emerges as the zero Kelvin extrapo- 
late from the PIMC simulatio ns as we have shown earlier 
(jKvlanpaa and Rantala|[20lol ). There, we evaluated the 
differences between three models for the description of 
the nuclear dynamics: Born-Oppenheimer approxima- 
tion, nuclei in thermal motion and nuclei with both ther- 
mal and quantum dynamics. At low temperatures the 
necessity of the fully quantum mechanical approach for 
all five particles was established. 

In the next section we present the essential details of 
the Feynman path integral quantum statistical approach, 
numerical simulation method and the model of the 
ion. In section 3 we present and analyze the energetics, 
partition function and other thermodynamic functions of 
the system using analytical forms where pertinent. In the 
last section the conclusions are given. 

2. METHOD AND MODELS 

According to the Feynman path in tegral formulati on of 
the quantum statistical mechanics ("Feynman '1998') the 
partition function of interacting distinguishable particles 

is given by the trace of the density matrix p{/3) = e~^^ 
as 

. AI-l 

Z=Tvp{(i)= URoARi. . . di?M-i n^^^^'^'''^'^'''^ 

where the action S{Ri, Ri+i]T) is taken over the path 
Ri — > Ri+i in imaginary time r — /3/M, where /? = 
l/fceT and M is called the Trotter number. The trace 
implies a closed path (Rm = Rq)- 

For s imulation we use the pair a pproximation in the 
action (|Storeiill968t lCeperlevlfT995l) for the Coulomb in- 
teraction of charges. This is exact in the limit M ^ oo, 
however, chemical accuracy is reached with sufficiently 
large M, i.e. small enough r. Sampling in the con- 
figuration space {Ri}^ in NVT ens emble is carried 
out using the Metropolis algorithm (jMetropolis et al\ 
119531) with bisection m oves and displacement moves 
(jChakravartv et aLl[T998l ). The total e nergy is calculated 
using the virial estimator ([Herman et al. 1982). 

The error estimate in the PIMC scheme is commonly 
given in power s of the imaginary time time-step r 
(|Ceperlevlll995[ ). Therefore, in order to systematically 
determine the thermal effects on the system we have car- 
ried out all the simulations with r ~ 0.03i?jj^, where E'h 
denotes the unit of Hartree. Thus, the temperatures and 
the Trotter number M are related by T = (fceMr)"^, 
where /cb is the Boltzmann constant. 

In the following we mainly use the atomic units, where 
the lengths, energies and masses are given in the units 
of Bohr radius (ao), Hartree (-Eh) and free electron mass 
(mo), respectively. Thus, for the mass of the electrons we 
take TTic = 1 and for the protons mp = 1.83615267248 x 
lO^TTio. Conversion of the units of energy is given by 
En = 219474.6313705cm~i « 27.2 eV, and correspond- 
ingly, fee = 3.1668152 x 10-^£;hK-^ 

The statistical standard error of the mean (SEM) with 
2 SEM limits is used as an error estimate for the evalu- 
ated observables. 

For the NVT simulations we place one Hjj^ ion, 
i.e. three protons and two electrons, into a cubic box and 
apply periodic boundary conditions and the minimum 



image principle. The simulations are performed in three 
different super cell (box) volumes: (300ao)3, (lOOao)^ 
and (50ao)^. These correspond to the mass densities 
of - 1.255 X 10^6 gcm-3, - 3.388 x 10"^ gcm-^ and 
~ 2.710 X 10"'* gcm~^, respectively. The density has no 
essential effect at low T, where dissociation rarely takes 
place. At higher T, however, the finite density gives rise 
to the molecular recombination balancing the more fre- 
quent dissociation. 

It should be pointed out that application of the mini- 
mum image principle with only one molecular ion in the 
periodic super cell may both give rise to the finite-size 
effects and also disregard higher density distribution ef- 
fects, i.e. fragments of several ions in the simulation box. 
Thus, the lower the density the better we minimize the 
finite-size effects, which in this work are negligible, if not 
absent. The zero density limit cannot be reached due to 
the finite T. To avoid all these ambiguities we define our 
targets as molecular energetics, molecular partition func- 
tion and other related molecular quantities. Therefore, 
in the following, we also exclude the trivial contribution 
from the center-of-mass thermal dynamics and energy 
^k-oT to the molecular quantities. 

The nuclear quantum dynamics, which was shown to 
be essential at low T, turns out to be negligible at higher 
temperatures. It is included, however, to be consistent 
with the low temperature results and our earlier study. 
For more details about the model and a discussion about 
the here n eglected contribution fr o m the exchange inter- 
action see IKvlanpaa and Rantalal (|2010f) . 

3. RESULTS AND DISCUSSION 

3.1. Overview of molecular energetics 

In Fig.[T]the NVT total energy (canonical ensemble in- 
ternal energy) of the H^ ion and its fragments is shown 
as a function of temperature. The molecular energy does 
not include the center-of-mass translational kinetic en- 
ergy ^k^T . The data from simulations are given as 
circles, squares and triangles corresponding to the three 
densities. The PIMC data is also given in Tables [T] and 

H 

The solid lines at T < 4000 K are fitted to analyti- 
cal model forms but at higher temperatures lines are for 
guiding the eye, only. Our low temperature fit and ana- 
lytical model, Eq. ([T]), is given as a blue dashed line and 
is discussed in the next sections in more detail. For com- 
pa rison the energies from the fitted partition function 
of iNeale and TennvsonI ()1995D is shown as black dots. 
These two do not manifest dissociation, and therefore, 
are not relevant at "higher T". 

The horizontal dash-dotted lines show the zero Kelvin 
energies for the ion and its fragments in Eq.([5]). One 
of these lines presents the energy for the "barrier to lin- 
earity", i.e. the energy needed for the transformation to 
linear molecular geometry. It is also seen to be roughly 
the barrier to dissociation within the considered molecu- 
lar densities. 

Above 4000 K the density dependence is clearly seen 
as varying composition of fragments. In the range from 
4000 to 10000 K the changing dissociation recombina- 
tion balance leads to distinctly different energetics, and 
above that, at our highest simulation temperatures the 
thermal ionization of hydrogen atoms starts contributing 
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Fig. 1. — (Color online) NVT total energy of the Hg molecu- 
lar ion as a function of temperature at three different densities: 
blue circles (~ 1.255 X 10~^ gcm""^), red squares (~ 3.388 X 
10~^ gcm~^) and gray triangles (~ 2.710 X 10~* gem"''). The 
blue dashed line is the energy fitted to Eq. (0. The black dots 
give the energy computed u sing the partition function fit given by 
INeale and TennvsonI II1995I) . The horizontal dash-dotted lines are 
the nonadiabatic zero Kelvin energies for the ion, its fragments 
and the barrier to linearity. The high temperature solid lines are 
mainly for guiding the eye, but used for numerical evaluation of 
the partition function, later. 

to the energy. However, it is worth pointing out that the 
temperature hmits of these three ranges, i.e. — 4000 
K, 4000 - 10000 K and above 10000 K, are subject to 
changes with larger variation of densities. 

Above 10000 K in our lowest density case the thermal 
ionization of H atoms is evident, see Fig. [1] but for our 
higher density cases some 15000 K is needed to show 
first signs o f ionization. Sirr ii lar tr end for the ionization 
is stated in iKoskinen et al\ (|2010D . although there the 
density is notably less than our lowest one. 

Let us now consider the dissociation recombination re- 
action chain, Eq. ([2]), and the contributing fragments 
to the quantum statistical NVT equilibrium trying to 
give an intuitive classical-like picture of the composition. 
With finite T, instead of zero, we have finite /3, instead 
of infinite, that brings classical nature to the system the 
more, the higher the temperature. In other words, the 
partial decoherence in our five particle quantum system 
increases with increasing temperature, that enables us 
to distinguish the fragments as separate molecules and 
atoms in thermal equilibrium. Based on this interpreta- 
tion, we show the total energy distribution in Fig. [5] from 
sampling the imaginary time paths at about 5000 K with 
M = 2048 in our highest simulation density. 

We see three main peaks and by inspection of the en- 
ergy distribution the first and the second can clearly be 
assigned to the rovibrationally excited and H2+H+, 
respectively. As expected, there are no rovibrational ex- 
citations available for 2H-t-H+, leading to the peak aver- 
age position very close to — IE'h- The fourth fragment. 




-1.3 -1.2 -1.1 -1.0 
E (units of Hartree) 



Fig. 2. — (Color online) Histogram of total energy sampling 
pinned in boxes of width O.OOIEh from (2 X 10*) X 10^ Monte Carlo 
samples averaged over blocks of 10^ samples. The simulation den- 
sity, temperature and Trotter number are ~ 2.710 X 10~* gcm~^, 
~ 5139.6 K and 2048, respectively. Other notations are taken from 

Fig.m 

H^+H, can be identified as the small high-energy side 
shoulder of H2+H"'" peak. With the interpretation of the 
area under the peak as the abundance of the fragment 
in the equilibrium we find this contribution to be much 
smaller than that of the others. 

It is important to note, however, that the above illus- 
tration is dependent on the block averaging procedure, 
see the caption of Fig. [21 Pinning the energy data of each 
and every sample, i.e. choosing block of size one sample, 
would broaden the peaks in Fig. [51 At the opposite limit, 
all samples in one block, the single mean energy or ensem- 
ble average is the quantum statistical expectation value 
-1.169(29)£;h, see Tabled Figs.[Iland[l where the sta- 
tistical uncertainty decreases with increasing simulation 
length. 

3.2. Molecular partition function 

To compare with the other published approaches for 
the molecular partition function based on single molecule 
quantum chemistry we start from the lowest temperature 
range from to ^ 4000 K, where the molecule does not 
essentially dissociate, yet. 

We present a low temperature Hj^ molecular partition 
function as a first approximation for the modeling of low 
density Hj^ ion containing atmospheres. Our aim is to 
find a simple analytical form, which can be accurately 
fitted to the NVT energies from our simulations. 

The partition function in terms of the Helmholtz free 
energy F is written as 

Z = e-P^, (3) 

where (3 = (fceT)"^, and the energy expectation value is 
straightforwardly derived from the partition function as 

<-> - -m 

After solving the free energy from Eq. ([3]) as 

F{T) ^-k^T In ZiT) (5) 



, : : , ^ 

2000 4000 6000 8000 10000 
T(K) 

Fig. 3.— (Color online) The molecular NVT ensemble InZ(T) 
from the energetics in Fig. [T] with the same notations. The blue 
solid line below 4000 K and its extrapolation (dashed line) are from 
Eq. JHJ, whereas the curves for three de nsities are from Eg. The 
black data (points) and fit (pluses) of iNeale and TennvsonrH1995l ) 
are also shown. 

we write F(T) — —k^Tf{T) and the energy expectation 
value may be written as 

(i?)=fcBT^^. (6) 

We find that a weU-behaving function fitting perfectly 
into our simulation data, 

{E) ^ knT^ {ae-''^ + c) + de-'^/^ , (7) 

allows analytical integration of Eq. ^ for /(T) or 
InZ(r), 

lnZ(r) = -|e-''^ + cT+-^e-"/^ + i?. (8) 

Using the boundary condition for the molecular partition 
function with a nondegenerate ground state, Z{Q) = 1 or 
InZ(O) — 0, we get D = a/h in our model. However, 
inclusion of the contributions from the ground state spin 
degeneracy factor and the zero-point rotations would give 
Z(0) = 1^ > 1, which would lead to D = a/h + In^, and 
thus, shift the function In Z by a constant, only. 

The weighted least squares fit of the above energy func- 
tion, Eq. ([7]), to our data for temperatures up to about 
3900 K, see Table [1] gives the parameters 

a = 0.00157426 
h = 0.000132273 
c = -6.15622 X 10"*^ 
d = 0.00157430 
a = 269.410 
D^a/hw 11.9016. 

In the fit, in addition to the (2SEM)^^ weights, we force 
the first derivative of the energy with respect to the tem- 
perature to be monotonically increasing up to 3900 K. 
The fit extrapolates the K energy to about 0.000549i?H 
above that of the para-H;^ , i.e. it gives an excellent match 
within the statistical error estimate. 

In Fig.[3]the function InZ(r) from Eq. ([5]) is shown in 
the range < T < 4000 K — the behaviour of the model 
at higher T is illustrated by the dashed line. Above 4000 
K the three curves for different densities are obtained 
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T(K) 

Fig. 4. — (Color online) Helmholtz free energy from Eq. JSJ (in 
units of Hartree). Notations are the same as in Fig. [3] 
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T(K) 

Fig. 5. — (Color online) Entropy from Eq. JTOj (in units of fee). 
Notations are the same as in Fig. [3] 

from those shown in Fig. [T] by numerical integration of 
Eq. dHJ as 

In Z{T) = In Z{T,) ^ d^, (9) 

wher e Ti = 500 K. 

Nea le and TennvsonI ()1995f ) have presented the parti- 
tion function In Z{T) based on a semi-empirical potential 
energy surface, see Fig. [3l The overall shape is similar 
to the one of ours. However, the energy (E) evaluated 
from their fit tends to be systematically lower than ours, 
although roughly within our 2 SEM error limits. Thus, 
the deviations are not visible in Fig. [T] For the parti- 
tion function the main difference is in their zero refer- 
ence, which also leads to In^ = — oo. This difference in 
zero reference goes back to that of energetics: our (E) 
in Eq. ^ a,t T — fits perfectly the zero Kelvin en- 
ergy of para- Hj", as mentioned earher. Th e zero Kelvin 
energy of the INeale and TennvsonI (jl995f) behind their 
partition function is that of the spin forbidden J — 
state, i.e. Z{0) — 0. This in fact, leads to the divergence 
of In Z{T) at T = 0. 

Our low temperature partition function Eq. ([5]) is close 
to complete. With the PIMC approach we implicitly in- 
clude all of the quantum states in the system with correct 
weight without any approximations. This partition func- 
tion is the best one for the modeling of the low density 
Hj^ ion containing atmospheres, at the moment. How- 
ever, it is valid up to about 4000 K, only. As soon as 
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1500 2000 2500 3000 3500 
T(K) 

Fig. 6. — (Color online) Heat capacity as a function of tempera- 
ture calculated using the analytical model of this work. The values 
on the y-axis are given in units of the Boltzmann constant k^. 

the density dependence starts playing larger role more 
complex models are needed. Such models can be fitted 
to our PIMC data given in Tables □ and [H 

3.3. Other thermodynamic functions 

In Fig. U we show the Helmholtz free energy from 
Eqs. ([5]) and dS]), combined. As expected, lower density 
or larger volume per molecule lowers the free energy due 
to the increasing entropic factor. Dissociation and the 
consequent fragments help in filling both the space and 
phase space more uniformly or in less localized manner. 

This kind of decreasing order is seen more clearly in 
the increasing entropy, shown in Fig. [S] The entropy has 
been evaluated from 

U -F 



S 



T 



(10) 



where the internal energy is U = (E) — {E)t=o- As 

expected, both the total energy (internal energy) and 
entropy reveal the dissociation similarly. 
Finally, in Fig. IH] we present the heat capacity 



Cv = 



d{E) 
dT ' 



(11) 



where (E) is taken from Eq. ([7]), which is valid at low 
temperatures, only. 

Considering the goodness of our functional form for 
(E), it is very convincing to see the plateau at about 
3/2 fee corresponding to "saturation" of the contribu- 
tion from the three rotational degrees of freedom. Thus, 
above 200 K the rotational degrees of freedom obey the 
classical equipartition principle of energy. It is the last 
term in the functional form of Eq. ([7]), that gives the 
flexibility for such detailed description of the energetics. 

It should be emphasized that the plateau is not ar- 
tificially constructed to appear at 3/2 fee, except for a 
restriction given for the first derivative of the total en- 
ergy to be increasing. Thus, the analytical model we 
present, Eq. ([7]), is found to be exceptionally successful 
at low temperatures. 

4. CONCLUSIONS 

We have evaluated the temperature dependent quan- 
tum statistics of the five particle molecular ion Hg far 



TABLE 1 

NVT ENERGETICS OF THE H'^ MOLECULAR, ION AT LOW 
TEMPERATURES — HERE THE SAME DATA APPLIES FOR ALL THREE 

DENSITIES. The energies are given in the units of Hartree 

(atomic units) and with 2 SEM ERROR ESTIMATES. THE 
ENERGIES FROM OUR LOW T FIT (LTFIT) FROM EQ. (0 AND 
THOSE FROM THE FIT Of INeALE AND TeNNYSOnI l|1995r (NT) ARE 
ALSO GIVEN AS COMPARISON. At K THE BEST UPPER BOUND IS 
GIVEN, SEE THE FOOTNOTE"^. 



TiK) 


PIMC' 


LTFIT'' 


NT fit'^ 







-1.3231 


(-1.32367)= 


~ 160.61 


-1.3227(7) 


-1.3227 


-1.3232 


~ 321.22 


-1.3221(6) 


-1.3220 


-1.3225 


~ 642.45 


-1.3198(6) 


-1.3202 


-1.3209 


~ 1052.6 


-1.3173(7) 


-1.3171 


-1.3179 


~ 1365.2 


-1.3143(5) 


-1.3141 


-1.3148 


~ 2000.3 


-1.3064(7) 


-1.3065 


-1.3070 


~ 2569.8 


-1.2983(8) 


-1.2984 


-1.2989 


~ 3049.2 


-1.2905(12) 


-1.2909 


-1.2917 


~ 3499.3 


-1.2840(12) 


-1.2835 


-1.2847 


~ 3855.6 


-1.2774(7) 


-1.2774 


-1.2792 



''This work. 

''Calculated from the fit given in lNeale and 'I'ennvsonI II1995I ). 
'^Para-Hg', see refs. IKutzelnieg and Jagueti II2006I) and 



IPavanello and Adamowiczl 120091) . 



TABLE 2 

PIMC NVT ENERGETICS OF THE 



MOLECULAR ION AT HIGH 



temperatures for the three densities (expressed as the 
number of molecular ions per volume), see fig. [t] 
Notations are the same as in Table [2] 



TiK) 


(300ao)"=* 


(100ao)-=* 


(50a„)--^ 


~ 3999.2 


-1.152(16) 






~ 4050.0 


-1.19(6) 






~ 4100.4 


-0.9995(4) 






~ 4498.2 


-0.9993(4) 


-1.219(34) 


-1.244(15) 


~ 4819.5 


-0.9993(4) 




-1.215(37) 


~ 5139.6 


-0.9995(4) 


-1.020(33) 


-1.169(29) 


~ 5634.8 






-1.156(66) 


~ 6070.3 


-0.9991(4) 


-1.018(18) 


-1.062(35) 


~ 7017.2 


-0.9995(4) 


-1.008(9) 


-1.024(12) 


~ 10279 


-0.997(3) 


-0.9995(8) 


-1.003(3) 


~ 12016 


-0.9993(6) 






~ 13997 


-0.86(10) 






~ 14951 


-0.805(23) 


-0.988(8) 


-0.9957(8) 



beyond its dissociation temperature at about 4000 K. 
This is done with the path integral Monte Carlo (PIMC) 
method, which is basis set and trial wave function free 
approach and includes the Coulomb interactions exactly. 
Thus, we are able to extend the traditional ab initio 
quantum chemistry with full account of correlations to 
finite temperatures without any approximations, also in- 
cluding the nuclear thermal and quantum dynamics. 

It is fair to admit, however, that PIMC is computa- 
tionally heavy for good statistical accuracy and approxi- 
mations are needed to solve the "Fermion sign problem" 
in cases where exchange interaction becomes essential. 

The temperature dependent mixed state description of 
the ion, the density dependent equilibrium dissocia- 
tion recombination balance and the energetics has been 
evaluated for the first time. With the rising temperature 
the rovibrational excitations contribute to the energet- 
ics, as expected, whereas the electronic part remains in 
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its ground state in the spirit of the Born-Oppenheimer 
approximation. At about 4000 K the fragments of the 
molecule, H2+H+, H^+H and 2H+H+, start contribut- 
ing. Therefore, Hg" ion becomes less dominant, and even- 
tually negligible in high enough T . 

We have shown how the partial dccohcrcncc in the 
mixed state can be used for interpretation of the frag- 
ment composition of the equilibrium reaction. Further- 
more, we have evaluated explicitly the related partition 
function, free energy, entropy and heat capacity, all as 
functions of temperature. An accurate analytical func- 



tional forms are given for the temperatures below disso- 
ciation. We consider all these as major improvements to 
the earlier published studies, where dissociation has not 
been considered. 
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